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The relation between force and stretch in the worm-like chain model of entropic elasticity is 
examined. Although no closed-form expression is valid for all values of forcing, solutions in the form 
of asymptotic series can be obtained under conditions of small and large applied force. The small and 
large stretch limits correspond to regular and boundary layer perturbation problems, respectively. 
The perturbation problems are solved and series solutions obtained for force as a function of stretch. 
The form of the asymptotic series suggest a uniform approximation valid for all stretch that is an 
improvement on existing approximations. 

PACS numbers: 87.15.-v, 46.15.Ff, 82.37.Rs, 87.16.Ac 
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I. INTRODUCTION 

The worm-like chain (WLC) is a model of entropic elas- 
ticity [1] for a macromolecule under thermal agitation. 
The main feature of the model, as compared to simpler 
ones such as the freely jointed chain (FJC) model Q, 
is the inclusion of bending energy. Applications of the 
WLC model range from macroscopic elasticity of rubber 
and elastomers Q to DNA unfolding Q|. With the in- 
crease in interest and application there is a need to more 
clearly understand how the WLC model relates mechan- 
ical parameters, and in particular, the relation between 
the force applied at the chain ends and the stretch. This 
is complicated by the implicit and complex functional 
dependence in the model. 

The objective of this paper is to provide, for the first 
time, explicit analytical expressions for the applied force 
as a function of the stretch of the WLC. We begin with 
a brief introduction of the WLC model, and a review of 
existing closed-form approximations to the force-stretch 
relationship. 



II. THE WORM-LIKE CHAIN MODEL 

An excellent overview of the theory underlying the 
WLC model is given by Marko and Siggia Q. Consider 
a uni-dimensional flexible chain of total length Lq with 
end-to-end applied force F. The free energy of the chain 
is 
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where L p is the persistence length, t(l) is the unit tangent 
vector, and (3 = (kT)^ 1 . The applied force results in 
average stretch z at temperature T. 

The natural non-dimensional units of force and stretch 



/ = PLpF, 



z/L . 



(2) 



Using standard arguments from statistical mechanics 0, 



L p dlnZ 



(3) 



where Z is the partition function over all possible states, 
it is certainly the case in elastomers, and generally true 
for DNA, that the persistence length is much less than the 
unfolded molecule end-to-end length. The large param- 
eter Lo/Lp ^> 1 ensures that InZ, which can be identi- 
fied as chain entropy, is dominated by the lowest energy 
state,. As a result @ Z —(L /L p )e a , where eo is a 
nondimcnsional energy, defined as 



eo = mm / da; 



1 
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The probability density function is normalized (ip, ip) 
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FIG. 1: The WLC relation between stretch s and applied force 
/. The numerical method is summarized in the Appendix. 
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with respect to the inner product 



(iM> 



dxip(x)(j)(x) 



(5) 



The function -0 is smooth and bounded for all — 1 < x < 
1. The stretch is then 



de 
df 



dxxtp . 



(6) 



The two terms in eo of (HJ) correspond to the bending and 
work terms in the original energy Ewlc, and the spe- 
cific form of the integrands is associated with rotational 
invariance about the force axis, with t-F = F cos 8 = Fx. 

The WLC problem therefore requires finding station- 
ary values of the functional 



ru-)= / dx±(l-x 2 W) 2 -f( J dxx^-s) 



e ( / dx^ 2 -l). 



(7) 



r(-0) contains the bending energy term plus two con- 
straints involving the first two moments of the function 
if). The normalization (tf>, ijj) = 1 defines -0 as a probabil- 
ity density function, while the constraint ^2 defines the 
stretch s. We may consider the stretch as given, so that / 
and eo are Lagrange multipliers, and the Euler-Lagrange 
equation is 

~ [(1 - x 2 W\' + fxiP + e V - 0, -1 < x < 1, (8) 

The objective is to find the lowest value of eo, and the 
force / is then uniquely determined as a function of s. 
This dictates an indirect procedure: consider / as given, 
and find eo, the lowest eigenvalue of the differential op- 
erator that depends upon /. Then s is determined as a 
function of / via either formulas given by eq. ©. Note 
that the value of T at the minimum is 70 = eo + s/, which 
is the Legendre transform of eo with / = d^fo/ds. The 
2D version of eq. ([5]) reduces to the Mathieu differential 
equation with solution in terms of Mathieu functions [7| . 
Prasad et al. [?J derived small and large force limits for 
the WLC in two dimensions using this approach. The 
focus here is on the 3D problem only. 

Figure [JJ shows the characteristic WLC curve, obtained 
from eqs. © and ((SJ) using a numerical method based 
on p|, see the Appendix. There are other ways to find 
/ = f(s), e.g. by solving the ODE using a shooting 
method [6]. The important issue is not, however, the 
numerical determination of the curve, but finding a suit- 
able analytic approximation. An excellent first step in 
this direction was made by Marko-Siggia [H who showed 
the leading order behavior for / <C 1 and for / 3> 1 
is / = |s and / _1 = 4(1 — s) 2 , respectively. Motivated 
by this limiting behavior they suggested the approximate 
functional form 
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This simple formula reproduces the small and large 
stretch leading order response in the respective limits. 
Ogden et al. || examined several alternative approxima- 
tions based on intelligent curve fitting to the f — s data 
in [6|. The simplest formula, which they called WLC3, is 
just the Marko-Siggia approximation with a single term 
added: 



WLCo = 



1 



1 



4(1 - s) 



9 71 s 7 s 

2 4 4 



(10) 



The extra quadratic term 



3 „2 



produces a dramatic im- 



provement, see Fig. [2j The root mean square error of 
WLC3 is 0.013 as compared with 0.339 for Jms- The an- 
alytical results of this paper will help explain this roughly 
25-fold increases in accuracy. We will return to consider 
WLC3 in Section |V] after deriving the small and large 
stretch approximations. The principal results of the pa- 
pers are summarized next. 



A. Summary of the main results 

The small and large stretch expansions are 

/ = 



3 „ 1 33 3 
2 * 20 S 



3393 5 
1400* 



1(1- s y> + m + mi 1 ~ s ) + 5i^( 1 ~ - s ) 2 + 



(11) 



valid for s <C 1 and 1 — s -C 1, respectively. Based on 
these limiting forms, and some numerical experimenta- 
tion, we find that the following approximation to / shows 
significant improvement on WLC3, 



WLCe = 



4^-i +S -^ 2 + ^ 3 (3-5 S )(19-20 S ). 

(12) 

This has rms error of 0.0047 and is compared with WLC3 
in Fig. H 

The remainder of the paper is organized as follows. 
The asymptotic series of eq. ([Til are derived in Sections 
Mil and IIVI The small stretch regime is considered first 
in Section IIII1 where the solution is developed using reg- 
ular perturbation methods. Large stretch is examined in 
Section IIVI Although the problem is a singular pertur- 
bation, it is reduced to a regular perturbation expansion 
using an inner scaled variable. The two asymptotic se- 
ries are compared with the exact solution in Section [V] 
The new and improved approximate formula valid for all 
values of stretch, large and small, is proposed after some 
numerical experimentation. 



III. SMALL STRETCH EXPANSION 

A. Perturbation theory 

Under small stretch, or equivalently small applied 
force, the WLC equation reduces to a regular pertur- 
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bation problem. Define 

L=^-(l-x^, (13) 
di dx 

then with the replacements eo —* 5 A and / — > the 
equation JHJ) becomes 

Lip + \ip + extp = 0, -l<x<l. (14) 

The small stretch limit corresponds to e -C 1. We seek 
solutions to eq. (114|) in the form of a regular perturbation 
expansion 

ip = ipo + eipi + e 2 ip 2 + • ■ ■ , (15a) 
A = A + eAi +e 2 A 2 + . . . . (15b) 

Substituting these into eq. (|14p and identifying terms 
of like order in the perturbation parameter e yields a 
sequence of equations. The first few of order e , e 1 and 
e 2 , are respectively, 





L ipo 


= 0, 


(16a) 




+ xipo + Xiipo 


= 0, 


(16b) 


L ip 2 + xtpi 4 


- AiV'i + \ 2 ipo 


= 0. 


(16c) 



where 

L = L + A . (17) 

Although the WLC corresponds to Ao = 0, it is useful 
to first consider the perturbation of an arbitrary ground 
state. 

The form of the 0(e fe ), k > 1, equation is 
L V>fc + xipk-i + Ai^ fe _i + \ 2 ipk-2 + ■ ■ ■ + Afe^o = 0- (18) 

The unperturbed solution tpo(x) is either an even or an 
odd function of x. It follows that ipk has the same or 



opposite parity depending as k is even or odd, respec- 
tively. We assume the unperturbed solution is normal- 
ized (ipo,ipo) = 1- 

The operator Lq is self adjoint with respect to the inner 
product ([5|), implying the solvability condition at 0(e fe ) 
is 

Afc + Afe_i(V>i,^o) + ■ ■ ■ + Ai(-0fc-i, ipo) + (xipk-i,ipo) = 0. 

The solvability condition essentially ensures that the so- 
lution to eq. (|18p can be expressed in terms of a sum of 
Legendre polynomials that are regular at the end points, 
i.e. P n . However, the expression for ipk has no compo- 
nent corresponding to ipo, in other words, (ipk, ipo) = $ko- 
Taking into account the parity of the successive terms 
gives the succinct result 

A a *_i=0, A afc = -<^fc_i,^d), fc = l,2,.... (19) 

Note that the first few equations simplify to 









L tpo 


= 0, 


(20a) 






Loipi 


+ xip 


= 0, 


(20b) 




L Q lp2 


4- xipi 4 


- A 2 ^o 


= 0, 


(20c) 




Lo4>z 


4- xip2 4 


- A 2 ^i 


= 0, 


(20d) 


L 1p4 - 


t- xip 3 4 


- A 2 -02 4 


- A 4 ^o 


= 0, 


(20e) 


Loips - 


1- X1p4 4 


- A 2 V'3 4 


- A 4 V'i 


= 0. 


(20f) 



We will solve these for the WLC problem, which corre- 
sponds to the lowest eigenvalue. Before considering the 
WLC specifically, we note some properties of the eigen- 
value perturbation that are valid for any eigenvalue. 

B. A2 for any initial state 

The unperturbed eigenvalue problem is Legendre's 
equation, and hence the most general form of the un- 
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perturbed solution is 

ipo(x) = c n P n (x), X = n(n + 1), (21) 
where P n is the Legendre polynomial of order n and the 
normalization factor is c n — 
Using the identity @ 

(2k + l)xP k = fcP fc _i + (fc + l)P fc+ i, (22) 
it is easy to show that 

L (P k+ i -P k -i(l-S k0 )) + 2(2fc + l)xP k = 0. (23) 
Hence, the first correction to the unperturbed mode is 

Cr, 



4>i 



2(2n + i; 



(Pn+l — -Pn-l(l — <^no))- 



(24) 



The first correction to the eigenvalue follows from the 
identities Q 



2n 



(xPj,P n ) = ^ggg ' 
t (2n+l)(2n+3) ' 



1, 



Z = 71 + 1, 



as 



A 2 = [2(2n- 1) (271 + 3)]- 



(25) 



(26) 



Note that A2 > for all n except n — 0, which has the 
lowest eigenvalue. We now consider the lowest energy 
state specifically and continue the perturbation expan- 
sion to higher orders. 



C. The lowest eigenvalue 

We focus on the unperturbed solution for n = 0, which 
has the lowest initial energy. The analysis of the previ- 
ous subsection gives the first two terms in the eigenvalue 
and eigenfunction expansions as Ao = 0, A2 = — i, and 

■0o = c oPo, "01 — QPii w ith Co = l/-\/2. These are the 
solutions of the first two in the hierarchy of equations 
(|2U)l . The next two are then solved to obtain 02 and 03, 
from which the next term in the eigenvalue expansion, 
A4, follows from eq. ([TO]) . 

In this manner the first six equations given in (|20[) may 
be solved successively. The terms in the eigenfunction 
expansion were obtained using Mathematica, 



ipo = co-Po, i>i = 5C0.P1, 02 = 7^2, 
z 18 



1p3 
"04 



CO 



3.4.5.6 

co 



7.8.9.10 

CO , 8 



(P 3 - HP), 



(aT-Ps 



212 p 



7520 



2 7 3 4 5 2 

and the corresponding expansion of the eigenvalue is 



(27a) 
(27b) 
(27c) 
(27d) 



A = -I £ 2 



11 



e 4 - 



47 



0(e 8 



(28) 



6 1080 34020 

The procedure can be continued; however the coefficients 
quickly become more unsightly. 



D. Small stretch expansion 

Taking into account the factor of 1/2 difference be- 
tween eq. (fl"4"|) and the WLC equation ([5]), the above 
analysis implies that the lowest perturbed energy is 



11 



3.47 



5.27 J 5.7.9.27 
The stretch follows from eq. ([6]), 



f 



3" 5.27 J 5.7.9.9 

and inverting the series gives 

, 3 33 o 9.13.29 - 

f = -s -\ s J H s° 

■' 2 20 1400 



(29) 



(30) 



(31) 



The accuracy of the small stretch expansion is shown 
in Fig. 02 with WLC3 used as a comparison. The relative 
error of the three term asymptotic series is less than 10 
for < s < 0.3, but the approximation deteriorates at 
higher values, as expected. 



IV. LARGE STRETCH: A BOUNDARY 
LAYER APPROXIMATION 

A. A singular perturbation problem 

The large stretch limit corresponds to large values of 
the applied force / in eq. ©. We therefore consider 



hLtp + Xip + e~ 2 x%jj = 0, -1<x<1, 



(32) 



for £ < 1. The second order differential operator L is 
defined in eq. (fTB"]) . and the factor of 1/2 is introduced 
for convenience. Equation (f3"2"| defines a singular per- 
turbation problem for tp(x), describing a boundary layer 
solution that is non-zero only near x = 1. In order to 
deduce this introduce the boundary layer variable 



X = (l-x)e~ 
Let = ip(x), and define A by 

A 



(33) 



(34) 



then the equation for ^ becomes 

{X^'Y + (A - X)^ - (e/2) (X 2 *')' = 0, (35) 

for < X < 2/e. This is now a regular perturbation 
problem in terms of the rescaled inner coordinate X. 
Note that the range of X depends upon the small pa- 
rameter e, although this is not a serious complication 
since the effective range of X is the positive real axis. 
Assuming the regular perturbation expansion 



* = *o + e*i +e 2 * 2 + 
A = A + eAi + e 2 A 2 + . 



(36a) 
(36b) 
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(a) (b) 

FIG. 3: The relative error of the perturbation solution compared with the exact solution at small stretch on a log scale (a) 
and absolute value (b). fs(s) is the three term expansion of eq. (|31[) . and /2(s) is the first two terms only. The relative error 
of the WLC3 approximation of eq. (llOf) is also shown. 



gives the sequence of equations 

{X%)' + (Ao - X)^ = 0, (37a) 
(X^Y + (Ao - X)*! + A^q - ±(A 2 * )' = 0, (37b) 

etc. The solution of the first equation, of order e°, is 

* (X) = C () e- X , A = 1, (38) 

where normalization implies Co = a/2- The next equa- 
tion, of order e , becomes 

(X^)' + (1 - X)*! + (A 1 +X- iX 2 )* - 0. (39) 

The solvability condition 




dXfAi + X - ±X 2 )* 2 (A) = 0, (40) 



implies the first correction is Ai = — j. 

It is evident that the solutions have the form of the fun- 
damental exponentially decaying solution ^q(X) multi- 
plied by polynomials in X. This suggests scaling ^ with 
respect to the leading order solution, 

*(X) =g{X)* {X). (41) 

The equation for g is 

Jg + eHg + (A-l - A 1 e)g = Q, (42) 

where the differential operators J and H are 

Jg(X) = Xg" + (1 - 2X)g\ (43a) 

Hg{X) = {X-^--\)g + (A 2 - X)g> ^-g". 

(43b) 



Assuming the expansion 

g = g + egi + e 2 g 2 + • • . , (44) 

then go = 1 and the equations for g\ through g± are 

J 9l + Hgo = 0, (45a) 

Jg 2 + Hgi + A 2 = 0, (45b) 

Jg 3 + Hg 2 + k 2 gi + A 3 = 0, (45c) 

J 54 + Hg 3 + A 2 g 2 + A 3gi + A 4 = 0, (45d) 

The procedure is then to find g\ as the particular solution 
to eq. (|45a[) and A 2 follows from the solvability condition 
for eq. (|45bj) : 




dX(A 2 + %)^(I) = 0. (46) 



These steps are repeated to find the successive functions 
gk and the eigenvalue coefficients A*,. 

Equations (|4"5|) were solved using Mathematica. We 
omit the detailed form of the gk functions and focus on 
the eigenvalue solution which is all that is required for 
the WLC model, 

1 1 1 1 3 2 885 3 „f 4<> tA* 

A = — -H e e 2 e 3 + 0(e 4 ). (47) 

e 2 e 4 64 512 262144 v 1 v ' 



B. Large stretch expansion 

The boundary layer solution with e = /~ 1//2 implies 
that the lowest energy state of the WLC has the large 
force expansion 

e »^ /+/2 l^-5W-ii4? + "- (48) 
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The stretch is then 
S = l 1 1 



5.9.59 



(49) 



2/3 2 7 /i 2 9 / 2 2 19 /^ 

Inverting the asymptotic series gives the desired expres- 
sion for / as a function of s, 

/= , 1 N9 +- + -(l- S ) + ^-(l-.s) 2 + .... (50) 
J 4(1- s) 2 32 64 v ' 32768 v ' v ' 

The large stretch asymptotic expansion is illustrated 
in Fig. [4j with WLC3 again used as a comparison. The 
relative error of the four term series is less than 10~ 3 for 
0.8 < s < 1, roughly. 



V. NUMERICAL EXPERIMENTS 

Comparison of the accuracy of the small and large 
stretch expansions in Figs. [3] and [4] indicate the at se- 
ries as developed here are accurate to within on part in 
10 3 for the range < s < 0.3 and 0.8 < s < 1, with 
zero error at s = and s = 1. This Section examines the 
question of finding an approximation that is uniformly 
valid over the entire range of stretch. 

The difference between the exact force function and 
WLC3 at small stretch follows from eqs. (fT0|) and (fTTj) 
as 



'26 
V 5 



10s 



1293 2 
175 A 



)+0(8% 



f - WLC 3 = 



j_ 

32 



(1-*) + 



27135 
32768 



(1 



(51) 



-0((1- 



The term —4s that distinguishes WLC3 from the 
Marko-Siggia approximation (5) therefore exactly can- 
cels the error in the latter at 0(s 2 ) in the small stretch 
limit. 



The quadratic for large stretch has zeros at s = 0.9191 
and s — 0.5337. The first zero, being close to s = 1 can 
be attributed as the cause of the zero of / — WLC3 at 
s i=s 0.9189, see Fig. (2). The zeros of the quadratic in 
(jSTj) for small stretch are complex. However, as Fig. (2) 
indicates / — WLC3 has a second zero at s w 0.5986. 
This property of WLC3, that it is exact at s ~ 0.6 and 
s « 0.92, partly explains its success as a uniform approx- 
imant. This suggests that any attempt at improving on 
WLC3 should maintain these zero crossings, and prefer- 
ably increase the number of zero crossings. 

At the same time wish to improve the accuracy at large 
stretch, requiring that the new approximation, say /*, is 
exact at s = 1. Consider the two parameter extension 
/* = WLC 3 +cs 3 {a-s), then the constraint /*(1) = 1/32 
implies c — (a — l)/32. Numerical experiments show 
that /* = WLC3 + |^ (a-i) ^ s n °t an improvement on 
WLC3 no matter what value of a is chosen. We therefore 
consider the two-parameter function 



/* = WLC 3 + 



s 3 (a — s)(b — s) 
32 (a - l)(b- 1)' 



(52) 



Using fminsearch in Matlab to minimize the root mean 
square error (/ — /*,/ — f*) 1 ^ 2 gives a = 0.5986, b — 
0.9458. Surprisingly, the value of a is precisely (to 
within four significant figures) the existing zero crossing 
of WLC3. In order to provide an approximation that is 
not too difficult to remember, we suggest rounding a and 
b up to 0.6 and 0.95, respectively. We call the resulting 
approximant WLCq, 

1 1 3 , 100 

WLC 6 = i(T -^- i+S - iS 2 + — S 3 (0.6- S )(0.95- S ). 

(53) 

The rms error incured by WLCq is 0.0047, as compared 
with 0.0045 for /* of (52) with a = 0.5986, b = 0.9458. 



10" 



of size (N + 1) x (JV + 1) with elements 



10" 



10" 



10"' 
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100 
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FIG. 5: The eigenvector components for / = 10 4 solved using 
Matlab with N = 300. Only the first 150 components are 
displayed, the rest are in the noise floor. 



The rms errors for /ms an d WLC3 are 0.3386 and 0.0132, 
respectively. These numbers indicate the remarkable ac- 
curacy of all three approximations to the exact force func- 
tion f(s). 



APPENDIX A: EXACT SOLUTION 

The exact solution for s = s(f) can be determined nu- 
merically quite easily [5[ . Define two symmetric matrices 



Sij = 



id i. 



V(2 i + l)(2j + l) : 



(Al) 



for i,j = 0, 1, 2, ... , N. Then for a given /, determine 
the minimum eigenvalue of D — /S and its eigenvector 
v. The strain is then 



s = v*Sv/(v'v). 



(A2) 



This algorithm can be effectively implemented in Matlab 
by using sparse matrix methods and the Matlab function 
eigs to find the single lowest eigenvalue. This is always 
negative but it is not always the smallest in magnitude, 
which is the criterion used in the function eigs. This can 
be circumvented by adding a multiple of the identity to 
D — /S so that the lowest eigenvalue is also the smallest 
in magnitude, without the eigenvector unchanged. We 
find that N=200 is more than sufficient to find s = s(f) 
f° r / < 10 with no apparent loss in numerical preci- 
sion. Figure [5] shows the amplitudes of the eigenvector 
components for / = 10 4 . Even at this large value the 
component with maximum amplitude is only v§. 
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